#include <string>
#include <algorithm>
#include <numeric>
#include <istream>
#include <vector>
#include <fstream>
#include <deque>
#include <cassert>
#include <map>
#include <set>
#include <stdexcept>
#include <ext/hash_map>
#include <ext/hash_set>


// boost functions
#include <boost/lexical_cast.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/lambda/lambda.hpp>

// my functions
#include <portable_timer.h>
#include <single_linkage.h>
#include <progress_indicator.h>
#include <blast_results.h>
#include <stlhelper.h>
#include <print_error.h>
#include <ultoa.h>
#include <index_string.h>
#include <single_linkage.h>

#include "get_arg.h"

#include <gc.h>
using namespace boost::lambda;
using std::accumulate;
using std::make_pair;
using std::map;
using std::cerr;
using std::ifstream;
using std::ofstream;
using std::pair;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::sort;
using boost::lexical_cast;
using std::ostringstream;


typedef __gnu_cxx::hash_map<unsigned, unsigned>	t_ks_index_per_gene;
typedef __gnu_cxx::hash_set<unsigned>			t_index_set;
typedef deque<t_ks_index_per_gene >				t_matrix_ks_indices;
typedef t_ks_index_per_gene::iterator			t_ks_index_iterator;
typedef t_ks_index_per_gene::const_iterator		t_ks_index_citerator;


//________________________________________________________________________________________

//		t_loci

//________________________________________________________________________________________
struct t_loci : public std::pair<unsigned, unsigned>
{
	t_loci():std::pair<unsigned, unsigned>(0U, 0U){}
	bool empty() const { return first == second;}
	bool assign(unsigned f, unsigned s)
	{
		bool was_empty = first == 0 && second == 0;
		first = f;
		second = s;
		return was_empty;
	}
};

bool loci_overlap(const t_loci& a, const t_loci& b)
{
	// if the end of one is completely beyond the opposite end of the other, no overlap
	return (a.empty() || b.empty() || !(a.second < b.first || b.second < a.first	));
}



//________________________________________________________________________________________

//		cluster_prot_by_loci_overlap

//			cluster transcripts if they overlap

//________________________________________________________________________________________

void cluster_prot_by_loci_overlap(	vector<unsigned>&			my_prot_ids,
									const vector<t_loci >& 		prot_id_to_loci,
									deque<vector<unsigned> >&	prot_clusters)
{
	typedef vector<unsigned>::iterator iter;
	//
	//	connect pairs which overlap
	//
	map<unsigned, unsigned> prot_ids_is_clustered;
	deque<pair<unsigned, unsigned> > overlapping_pairs;
	for (iter ii = my_prot_ids.begin(); ii != my_prot_ids.end(); ++ii)
		for (iter jj = my_prot_ids.begin(); ii != jj; ++jj)
		{
			assert(*ii < prot_id_to_loci.size());
			assert(*jj < prot_id_to_loci.size());
			if (loci_overlap(prot_id_to_loci[*ii], prot_id_to_loci[*jj]))
			{
				overlapping_pairs.push_back(make_pair(*ii, *jj));
				++prot_ids_is_clustered[*ii];
				++prot_ids_is_clustered[*jj];
			}
		}

	//
	//	single linkage cluster
	//
	if (overlapping_pairs.size())
	{
		single_linkage_clusters(overlapping_pairs, prot_clusters);
		assert(prot_clusters.size());
		assert(prot_clusters.size() <= my_prot_ids.size());
	}


	//
	// add clusters for each prot_id with no overlap at all
	//
	for (iter ii = my_prot_ids.begin(); ii != my_prot_ids.end(); ++ii)
		if (!prot_ids_is_clustered[*ii])
		{
			prot_clusters.push_back(vector<unsigned>());
			prot_clusters.back().push_back(*ii);
		}

	//
	// check all protein IDs in one cluster or another
	//
	#ifndef NDEBUG
	map<unsigned, unsigned> prot_ids_cnts_in_cluster;
	for (unsigned ii = 0; ii < prot_clusters.size(); ++ii)
		for (unsigned jj = 0; jj < prot_clusters[ii].size(); ++jj)
			++prot_ids_cnts_in_cluster[prot_clusters[ii][jj]];
	assert(prot_ids_cnts_in_cluster.size() == my_prot_ids.size());
	typedef map<unsigned, unsigned>::iterator miter;
	unsigned total = 0;
	for (miter i = prot_ids_cnts_in_cluster.begin(), e = prot_ids_cnts_in_cluster.end();
																			i != e; ++i)
		total += i->second;
	assert(total == my_prot_ids.size());
	#endif
}

//________________________________________________________________________________________

//		set_best_transcript_and_split_gene

//________________________________________________________________________________________
struct op_no_loci_overlap
{
	const t_loci& 					loci;
	const vector<t_loci >& 			prot_to_loci;
	op_no_loci_overlap(const t_loci& l_, const vector<t_loci >& prot_to_loci_)
							:loci(l_), prot_to_loci(prot_to_loci_){}
	bool operator()(unsigned prot_id)
	{	return !loci_overlap(loci, prot_to_loci[prot_id]);}
};

void set_best_transcript_and_split_gene(unsigned						gene_id,
										unsigned						best_prot_id,
										deque< vector<unsigned> >&		gene_id_to_prot_id,
										t_index_string&					gene_names,
										const vector<t_loci >& 			prot_id_to_loci,
										map<unsigned, unsigned>&		cnt_split_genes,
										vector<unsigned>&				extra_genes)
{
	vector<unsigned> prot_ids;
	typedef vector<unsigned>::const_iterator citer;
	citer i		= gene_id_to_prot_id[gene_id].begin();
	citer end	= gene_id_to_prot_id[gene_id].end();
	copy_if(gene_id_to_prot_id[gene_id], prot_ids,
				op_no_loci_overlap(prot_id_to_loci[best_prot_id], prot_id_to_loci));

	// check if divide into one or more
	deque<vector<unsigned> >	prot_clusters;
	cluster_prot_by_loci_overlap(prot_ids, prot_id_to_loci, prot_clusters);
	assert(prot_clusters.size());
	if (prot_clusters.size() != 1)
	{
		// add names of split gene
		unsigned base_gene_id = gene_names.size();
		string gene_name = gene_names.string_at(gene_id);
		for (unsigned i = 0; i < prot_clusters.size(); ++i)
		{
			unsigned index= ++cnt_split_genes[gene_id];
			extra_genes.push_back(gene_names.index(gene_name + "_EXTRA_" + ultoa(index)));
			gene_id_to_prot_id.push_back(vector<unsigned>());
			gene_id_to_prot_id.back().swap(prot_clusters[i]);
		}

		#ifndef NDEBUG
		cout << gene_names.string_at(gene_id) << "\n";
		#endif

	}

	gene_id_to_prot_id[gene_id].clear();
	gene_id_to_prot_id[gene_id].push_back(best_prot_id);
}



//________________________________________________________________________________________

//		open_existing_file

//________________________________________________________________________________________
unsigned open_existing_file(ifstream& is, const string& file_name)
{
	namespace fs = boost::filesystem;
	if (!fs::exists(file_name))
		throw std::runtime_error(file_name + " does not exist.");
	if (fs::is_directory(file_name))
		throw std::runtime_error(file_name + " is a directory not a file.");
	boost::intmax_t file_size = fs::file_size( file_name);
	if (!is)
		throw std::runtime_error("Could not open " + file_name);
	return file_size;
}


//________________________________________________________________________________________

//		lookup_prot_name

//________________________________________________________________________________________
inline unsigned lookup_prot_name(const t_index_string& prot_names, const string& prot_name)
{
	unsigned prot_id;
	if (!prot_names.is_indexed(prot_name, prot_id))
		throw std::runtime_error("Protein identifier in ks result [" + prot_name +
								"] not found in the supplied protein id to gene id "
								"lookup table.");
	return prot_id;
}





//________________________________________________________________________________________

//		parse_ks_res_from_file

//________________________________________________________________________________________
void parse_ks_res_from_file(	const string&			file_name,
								const deque<unsigned>&	prot_id_to_gene_id,
								const t_index_string&	prot_names,
								const unsigned			cnt_genes_ids,
								const deque<char>&		single_transcript_gene_ids,

								// multiple transcript genes saved
								deque< vector<unsigned > >&	gene_id_to_prot_id,
								deque<float>&			ks_res,
								t_matrix_ks_indices& 	matrix_ks_indices,

								t_progress_indicator& 	dots,
								t_program_args&			args)
{

	unsigned				cnt_multi_transcript_gene_results = 0;
	vector<char>			gene_ids_with_ks(cnt_genes_ids, 0);
	vector<char>			prot_ids_with_ks(prot_id_to_gene_id.size(), 0);


	const unsigned DOTS_SIZE = 20000;
	dots.set_every(DOTS_SIZE);

	unsigned cnt_lines = 0;

	ifstream is;
	unsigned total_file_size = open_existing_file(is,  file_name);

	try
	{

		while (1)
		{
			// get data for one protein pair
			string prot_name1, prot_name2;
			is >> prot_name1 >> prot_name2;
			string ignore_ka;
			is >> ignore_ka;
			float ks;
			is >> ks;
			is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
			if (!is.good())
				break;

			if (dots.value() % (DOTS_SIZE * 10) == 0)
			{
				boost::intmax_t file_size_so_far = is.tellg();
				unsigned dots_needed = unsigned(dots.value() * 1.0 *
												total_file_size / file_size_so_far);
				dots_needed += DOTS_SIZE;
				dots.set_limit(dots_needed);
			}
			cnt_lines++;
			dots++;

			// get protein and gene indices
			unsigned prot_id1 = lookup_prot_name(prot_names, prot_name1);
			unsigned prot_id2 = lookup_prot_name(prot_names, prot_name2);
			unsigned gene_id1 = prot_id_to_gene_id[prot_id1];
			unsigned gene_id2 = prot_id_to_gene_id[prot_id2];

			// note the gene_names
			gene_ids_with_ks[gene_id1] = 1;
			gene_ids_with_ks[gene_id2] = 1;

			// single transcript gene only have one choice of transcript: ignore
			if (single_transcript_gene_ids[gene_id1] &&
				single_transcript_gene_ids[gene_id2])
				continue;

			// if first time encountered this protein, add to list of proteins in gene
			if (!prot_ids_with_ks[prot_id1])
			{
				gene_id_to_prot_id[gene_id1].push_back(prot_id1);
				prot_ids_with_ks[prot_id1] = 1;
			}
			if (!prot_ids_with_ks[prot_id2])
			{
				gene_id_to_prot_id[gene_id2].push_back(prot_id2);
				prot_ids_with_ks[prot_id2] = 1;
			}

			// save ks results and its index for prot1 vs prot2 if not already saved
			if (matrix_ks_indices[prot_id1].count(prot_id2) == 0)
			{
				matrix_ks_indices[prot_id1][prot_id2] =
				matrix_ks_indices[prot_id2][prot_id1] = ks_res.size();
				ks_res.push_back(ks * ks * ks);
				++cnt_multi_transcript_gene_results;
			}
		}

		if (is.eof())
		{
			unsigned cnt_gene_ids_with_results = accumulate(gene_ids_with_ks.begin(),
																gene_ids_with_ks.end(),  0U);
			args.VERBOSE() << "\t"	<< cnt_gene_ids_with_results
									<< " genes with dS results in total\n";
			args.VERBOSE() << "\t"	<< dots.value() << " dS results in total\n";
			args.VERBOSE() << "\t"	<< cnt_multi_transcript_gene_results
									<< " dS results for genes with multiple transcripts.\n";

			return;
		}
	}
	catch(std::runtime_error& e)
	{
		throw std::runtime_error("Invalid ks results format on line #" +
									ultoa(cnt_lines) + " of " + file_name + "\n" + e.what());
	}
	throw std::runtime_error("Invalid ks results format on line #" +
								ultoa(cnt_lines) + " of " + file_name);


}


//________________________________________________________________________________________

//		parse_protid_loci_from_file

//________________________________________________________________________________________
void parse_protid_loci_from_file(	t_index_string&			prot_names,
									vector<t_loci>&			prot_id_to_loci,
									t_progress_indicator& 	dots,
									t_program_args&			args)
{
	ifstream is;
	open_existing_file(is, args.f_prot_loci);

	unsigned cnt_lines = 0;
	unsigned cnt_prot_loci = 0;
	unsigned cnt_prot_loci_ignored = 0;
	unsigned cnt_prot_loci_equal = 0;

	// get data for one protein
	string prot_name;
	unsigned start, finish;
	dots.set_every(20000);
	while (1)
	{
		// get next line
		is >> prot_name >> start >> finish;
		is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
		if (!is.good())
			break;

		dots++;
		cnt_lines++;
		unsigned prot_id;
		if (!prot_names.is_indexed(prot_name, prot_id))
		{
			++cnt_prot_loci_ignored;
			continue;
		}
		if (start == finish)
		{
			++cnt_prot_loci_equal;
			continue;
		}

		// make sure right way around
		if (start > finish)
			std::swap(start, finish);

		// save
		assert(prot_id < prot_id_to_loci.size());
		if (prot_id_to_loci[prot_id].assign(start, finish))
			++cnt_prot_loci;

	}

	if (is.eof())
	{
		dots.finish();
		args.VERBOSE() << "\tFinished parsing\t" << cnt_prot_loci << " transcript loci\n";
		if (cnt_prot_loci_ignored)
			args.VERBOSE() << "\t" << cnt_prot_loci_ignored << " transcript loci had no "
										"gene identifier mapping and were discarded.\n";
		if (cnt_prot_loci_equal)
			args.VERBOSE() << "\t" << cnt_prot_loci_ignored << " transcripts had empty loci "
										" (start = end) and were discarded.\n";
		return;
	}
	throw std::runtime_error("Invalid protein name and loci on line #" +
								ultoa(cnt_lines) + " of " + args.f_prot_loci);


}


//________________________________________________________________________________________

//		parse_protid_gene_id_from_file

//________________________________________________________________________________________
void parse_protid_gene_id_from_file(const string&		  file_name,
									deque<unsigned>&	  prot_id_to_gene_id,
									deque<unsigned>&	  prot_id_to_species_id,
									t_index_string&		  prot_names,
									t_index_string&		  species_names,
									t_index_string&		  gene_names,
									deque<char>&		  single_transcript_gene_ids,
									t_progress_indicator& dots,
									t_program_args&		  args)
{
	ifstream is;
	open_existing_file(is, file_name);

	unsigned cnt_lines = 0;

	dots.set_every(20000);
	while (1)
	{
		// get data for one gene
		string prot_name, gene_name, species_name;
		is >> prot_name >> gene_name >> species_name;
		is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
		if (!is.good())
			break;
		dots++;
		cnt_lines++;

		bool added;
		unsigned prot_id = prot_names.index(prot_name, added);
		unsigned species_id = species_names.index(species_name);
		if (added)
		{
			bool new_gene_id;
			unsigned gene_id = gene_names.index(gene_name, new_gene_id);

			// has a transcript for this gene been encountered before.
			// if not, it may be a single transcript gene
			if (new_gene_id)
			{
				assert(gene_id == single_transcript_gene_ids.size());
				single_transcript_gene_ids.push_back(1);
			}
			else
				// 2nd transcript encountered: must be multiply spliced
				single_transcript_gene_ids[gene_id] = 0;

			// save prot_id lookup
			prot_id_to_gene_id.push_back(gene_id);
			prot_id_to_species_id.push_back(species_id);

		}

		// prot_id shouldn't appear twice in this table
		else
		{
			args.CERR() << "WARNING:\tDuplicated prot ID " << prot_name
						<< " with previous gene of " << gene_names.string_at(prot_id_to_gene_id[prot_id])
						<< " and now with gene of " << gene_name << "\n";
		}
	}

	if (is.eof())
	{
		dots.finish();
		args.VERBOSE() << "\tFinished parsing all possible identifiers:\n";
		args.VERBOSE() << "\t\t"	<< prot_names.size()
									<< " possible protein transcript identifiers from \n";
		args.VERBOSE() << "\t\t"	<< species_names.size() << " species.\n";
		args.VERBOSE() << "\t\t"	<< gene_names.size()
									<< " possible gene identifiers, of these \n";
		assert(gene_names.size() == single_transcript_gene_ids.size());
		unsigned cnt_single_transcript_gene =
			accumulate(single_transcript_gene_ids.begin(),  single_transcript_gene_ids.end(), 0U);
		args.VERBOSE() << "\t\t"	<< cnt_single_transcript_gene
									<< " are for genes with single transcripts.\n";
		return;
	}

	throw std::runtime_error("Invalid protein / gene name on line #" +
								ultoa(cnt_lines) + " of " + file_name);


}




//________________________________________________________________________________________

//		link_genes_by_closest_ortholog

//________________________________________________________________________________________
void link_genes_by_closest_ortholog(const deque<float>&		   	ks_res,
									const t_matrix_ks_indices&	matrix_ks_indices,
									const deque<unsigned>&		prot_id_to_species_id,
									const deque<unsigned>&		prot_id_to_gene_id,
									const unsigned				cnt_species,
									deque<vector<unsigned> > &	gene_clusters
									)
{
	deque<pair<unsigned, unsigned> > neighbouring_gene_pairs;
	const unsigned undef = std::numeric_limits<unsigned>::max();
	vector<unsigned> neighbouring_prot_ids;
	vector<float> neighbouring_ks;
	for (unsigned prot_id = 0; prot_id < matrix_ks_indices.size(); ++prot_id)
	{
		// each species has a closest neighbour to this
		neighbouring_prot_ids.assign(cnt_species, undef);
		neighbouring_ks.assign(cnt_species, 10000);

		t_ks_index_citerator i = matrix_ks_indices[prot_id].begin();
		t_ks_index_citerator end = matrix_ks_indices[prot_id].end();

		for (; i != end; ++i)
		{
			unsigned other_prot_id = i->first;
			unsigned other_species_id = prot_id_to_species_id[other_prot_id];
			float other_ks = ks_res[i->second];
			assert(other_species_id < cnt_species);

			// save if smallest Ks for this species
			if (other_ks < neighbouring_ks[other_species_id])
			{
				neighbouring_ks[other_species_id] = other_ks;
				neighbouring_prot_ids[other_species_id] = prot_id;
			}
		}
		// ignore same species comparisons (want closest orthologues only)
		neighbouring_prot_ids[prot_id_to_species_id[prot_id]] = undef;

		// if there is a neighbour, look up its gene and add gene pair to neighbour
		for (unsigned i = 0; i < neighbouring_prot_ids.size(); ++i)
			if (neighbouring_prot_ids[i] != undef)
				neighbouring_gene_pairs.push_back(
						std::make_pair(	prot_id_to_gene_id[prot_id],
										prot_id_to_gene_id[neighbouring_prot_ids[i]]));

	}

	single_linkage_clusters(neighbouring_gene_pairs, gene_clusters);

}

//________________________________________________________________________________________

//		resolve_best_transcripts

//________________________________________________________________________________________
float get_ks(	unsigned					i,
				unsigned					j,
				const deque<float>&			ks_res,
				const t_matrix_ks_indices&	matrix_ks_indices)
{
	if (matrix_ks_indices[i].count(j) == 0)
		return 125.0F;
	else
		// convoluted syntax because array ref operator can insert, but using const object
		return ks_res[matrix_ks_indices[i].find(j)->second];
}
struct t_recurse_find_best_transcript
{
	float								best_ks_total;
	vector<unsigned>					best_multi_prot_ids;
	vector<unsigned>					singles_prot_ids;
	vector<unsigned>					multi_gene_ids;
	const deque<float>&					ks_res;
	const t_matrix_ks_indices&			matrix_ks_indices;
	const deque< vector<unsigned> >&	gene_id_to_prot_id;
	//
	//	reset for
	//
	void								next_round()
	{
		best_ks_total = std::numeric_limits<float>::max();
		best_multi_prot_ids.clear();
		multi_gene_ids.clear();
		singles_prot_ids.insert(singles_prot_ids.end(),
								best_multi_prot_ids.begin(),
								best_multi_prot_ids.end());
	}
	t_recurse_find_best_transcript(	const deque<float>&				ks_res_,
									const t_matrix_ks_indices&		matrix_ks_indices_,
									deque< vector<unsigned> >&		gene_id_to_prot_id_)
								:ks_res(ks_res_), matrix_ks_indices(matrix_ks_indices_),
									gene_id_to_prot_id(gene_id_to_prot_id_),
									best_ks_total(std::numeric_limits<float>::max()){}

	void recurse(unsigned level, vector<unsigned>& multi_prot_ids, float ks);
};

void t_recurse_find_best_transcript::recurse(unsigned level,
											vector<unsigned>& multi_prot_ids, float ks)
{
	if (level > multi_gene_ids.size())
	{
		if (ks < best_ks_total)
		{
			best_ks_total = ks;
			best_multi_prot_ids = multi_prot_ids;
			return;
		}
	}

	unsigned gene_id = multi_gene_ids[level];
	for (unsigned i = 0; i < gene_id_to_prot_id[gene_id].size(); ++i)
	{
		unsigned prot_id = gene_id_to_prot_id[gene_id][i];
		float curr_ks = ks;
		for (unsigned j = 0; j < singles_prot_ids.size(); ++j)
			 curr_ks += get_ks(singles_prot_ids[j], prot_id, ks_res, matrix_ks_indices);
		for (unsigned j = 0; j < multi_prot_ids.size(); ++j)
			 curr_ks += get_ks(multi_prot_ids[j], prot_id, ks_res, matrix_ks_indices);
		multi_prot_ids.push_back(prot_id);
		recurse(level + 1,  multi_prot_ids,  curr_ks);
		multi_prot_ids.pop_back();
	}
}

void resolve_best_transcripts(	vector<unsigned>&				gene_cluster,
								const deque<float>&				ks_res,
								const t_matrix_ks_indices&		matrix_ks_indices,
								const vector<t_loci >& 			prot_id_to_loci,
								deque< vector<unsigned> >&		gene_id_to_prot_id,
								map<unsigned, unsigned>&		cnt_split_genes,
								t_index_string&					gene_names)
{
	//
	//	divide into multiple and single transcript genes
	//
	vector<unsigned> gene_singles;
	t_recurse_find_best_transcript calc(ks_res, matrix_ks_indices, gene_id_to_prot_id);


	typedef vector<unsigned>::const_iterator citer;
	for (citer gi = gene_cluster.begin(), end = gene_cluster.end(); gi != end; ++gi)
	{
		if (gene_id_to_prot_id[*gi].size() == 1)
		{
			gene_singles.push_back(*gi);
			calc.singles_prot_ids.push_back(gene_id_to_prot_id[*gi][0]);
		}
		else
			calc.multi_gene_ids.push_back(*gi);
	}

	while (1)
	{
		//
		//	split genes into singles and multis
		//
		typedef vector<unsigned>::const_iterator citer;
		for (citer gi = gene_cluster.begin(), end = gene_cluster.end(); gi != end; ++gi)
		{
			if (gene_id_to_prot_id[*gi].size() == 1)
			{
				gene_singles.push_back(*gi);
				calc.singles_prot_ids.push_back(gene_id_to_prot_id[*gi][0]);
			}
			else
				calc.multi_gene_ids.push_back(*gi);
		}

		//
		//	if no more multis, finish
		//
		if (calc.multi_gene_ids.empty())
			break;

		//
		// recursively find best transcripts
		//
		vector<unsigned> multi_prot_ids;
		calc.recurse(0, multi_prot_ids, 0.0f);

		//
		// assign best transcripts per gene and split gene if necessary
		//
		vector<unsigned> extra_genes;
		for (unsigned i = 0; i < calc.multi_gene_ids.size(); ++i)
			set_best_transcript_and_split_gene(	calc.multi_gene_ids[i],
												calc.best_multi_prot_ids[i],
												gene_id_to_prot_id, gene_names,
												prot_id_to_loci, cnt_split_genes,
												extra_genes);

		//
		// all previouly multi genes have been resolved and are now singles
		//
		gene_singles.insert(	gene_singles.end(),
								calc.multi_gene_ids.begin(), calc.multi_gene_ids.end());
		calc.multi_gene_ids.clear();
		calc.next_round();

		//
		//	check multis and singles for extra genes
		//
		gene_cluster.swap(extra_genes);
	}
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	// get command line arguments
	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;

	//
	//	protein to gene id lookup
	//
	deque<unsigned>						prot_id_to_gene_id;
	t_index_string						prot_names;
	t_index_string						gene_names;
	t_index_string						species_names;
	deque<char>							single_transcript_gene_ids;
	deque<unsigned	>					prot_id_to_species_id;
	{
		args.VERBOSE() << "\tReading gene to protein mapping...\n";
		t_progress_indicator 				dots(args.VERBOSE(),  5000);
		parse_protid_gene_id_from_file(	args.f_prot_to_gene_id,
										prot_id_to_gene_id,
										prot_id_to_species_id,
										prot_names,
										species_names,
										gene_names,
										single_transcript_gene_ids,
										dots,
										args);
	}


	//
	//	protein to genomic loci lookup
	//
	vector<t_loci> prot_id_to_loci(prot_names.size());
	{
		args.VERBOSE() << "\tReading protein loci...\n";
		t_progress_indicator 				dots(args.VERBOSE(),  5000);

		parse_protid_loci_from_file(prot_names, prot_id_to_loci, dots, args);
		dots.finish();
	}

	deque< vector<unsigned> >	gene_id_to_prot_id(gene_names.size());
	deque<float>				ks_res;
	t_matrix_ks_indices			matrix_ks_indices;
	{
		args.VERBOSE() << "\tParsing blast results from [" << args.f_ks << "]. \n";
		t_progress_indicator 				dots(args.VERBOSE(),  5000);
		dots.use_timer();

		parse_ks_res_from_file(	args.f_ks,
								prot_id_to_gene_id,
								prot_names,
								gene_names.size(),
								single_transcript_gene_ids,

								// multiple transcript genes saved
								gene_id_to_prot_id,
								ks_res,
								matrix_ks_indices,
								dots,
								args);
	}

	// sort prot_ids (smaller protein ids always preferred);
	// 	i.e. if you sort prot to gene_id table in reverse transcript size, then large
	//	transcripts preferred
	for (unsigned i = 0; i < gene_id_to_prot_id.size(); ++i)
		sort(gene_id_to_prot_id[i].begin(), gene_id_to_prot_id[i].end());

	//
	//	Find the closest neighbour (ks) in the other species for each transcript
	//	and single linkage cluster genes
	deque<vector<unsigned> > gene_clusters;
	link_genes_by_closest_ortholog(ks_res,
								   matrix_ks_indices,
								   prot_id_to_species_id,
								   prot_id_to_gene_id,
								   species_names.size(),
								   gene_clusters);



	//
	// for each gene cluster resolve best transcript
	//
	map<unsigned, unsigned>		cnt_split_genes;
	#ifndef NDEBUG
	unsigned orig_gene_count = gene_names.size();
	#endif
	for (unsigned i = 0; i < gene_clusters.size(); ++i)
		resolve_best_transcripts(gene_clusters[i], ks_res, matrix_ks_indices, prot_id_to_loci,
									gene_id_to_prot_id, cnt_split_genes, gene_names);
	assert(cnt_split_genes.size() == gene_names.size() - orig_gene_count);

	//
	//	print out gene to prot ids
	//

	assert(gene_names.size() == gene_id_to_prot_id.size());
	for (unsigned i = 0; i < gene_id_to_prot_id.size(); ++i)
	{
		if (gene_id_to_prot_id[i].size() == 1)
		{
			assert(gene_id_to_prot_id[i].size() == 1);
			cout	<< gene_names.string_at(i)	<< "\t"
					<< prot_names.string_at(gene_id_to_prot_id[i][0])	<< "\n";
		}
	}


	args.VERBOSE() << "\tAll finished after " << timer.pretty_reset() << "\n";
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}



